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ABSTRACT 

The interaction of optically emitting clouds with warm X-ray gas and hot, tenuous 
radio plasma in radio jet cocoons is modelled by 2D compressible hydrodynamic sim- 
ulations. The initial setup is the Kelvin-Helmholtz instability at a contact surface of 
density contrast 10 4 . The denser medium contains clouds of higher density. Optically 
thin radiation is realised via a cooling source term. The cool phase effectively extracts 
energy from the other gas which is both, radiated away and used for acceleration of 
the cold phase. This increases the system's cooling rate substantially and leads to a 
massively amplified cold mass dropout. We show that it is feasible, given small seed 
clouds of order 100 Mq, that all of the optically emitting gas in a radio jet cocoon may 
be produced by this mechanism on the propagation timescale of the jet. The mass is 
generally distributed as T -1 / 2 with temperature, with a prominent peak at 14,000 K. 
This peak is likely to be related to the counteracting effects of shock heating and a 
strong rise in the cooling function. The volume filling factor of cold gas in this peak 
is of the order 10 -5 to 10 -3 and generally increases during the simulation time. 
The simulations tend towards an isotropic scale free Kolmogorov-type energy spec- 
tr um over the simulation tim escale. We find the same Mach- number density relation 
as iKritsuk fc Normanl (|2004l ) and show that this relation may explain the velocity 
widths of emission lines associated with high redshift radio galaxies, if the environ- 
mental temperature is lower, or the jet-ambient density ratio is less extreme than in 
their low redshift counterparts. 

Key words: hydrodynamics - instabilities - turbulence - galaxies: jets - methods: 
numerical. 



1 INTRODUCTION 

Extragalactic radio jets have been observed to interact 
strongly with their environment. This is evident not only 
from their often fat radio cocoons but also from the cospatial 
optical and X-ray emission (detailed below). This is observa- 
tional evidence for what is expected from theory, namely the 
co-existence of hot magnetised cocoon plasma, warm X-ray 
gas, and cold, possibly star forming clouds, all interacting 
with one another. 

This article is particularly stimulated by obser- 
vations of radio galaxies at a redshift of about one (e.g. 
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served: the smaller radio sources have the larger emission 
line regions, the greater line widths, and are more consistent 
with shock ionisation. This is further evidence for the tight 
connection between the gas phases present. 

The traditional viewpoint has been to study the indi- 
vidual phases separately, analog ously to interstellar medium 
(ISM) turbulence theory (|Eimegreen fc Scald l2004h . As for 
the case of the ISM, new insights can be expected from direct 
hydrodynamical simulation of the interaction. These inter- 
actions are so far unknown. The open questions include en- 
ergy, momentum, and mass transfer between the phases. The 
origin of the c old phase is unclear (apart from the merger 
cases, compare INeeser et al.lll999T ): a common idea is that 
pre-existing clouds are activated for example by the bow 
shock or photo-ionised by a central quasar. A possible prob- 
lem of this model might be the low visibility of clouds far 
from the radio e mission, which wou ld still be expected to 
be photo-ionised (|lnskip et al.l l2002bh . It is rather unlikely 
that much of the cool gas is drawn from the host galaxy, 
because the solid angle occupied by the beam is quite small, 
and except at the very end of the jet the densest clouds are 
likely to pass through the contact surface into the cocoon, 
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where the fl ow if anything moves the c louds back towards the 
galaxy fe.g. I Alexander &; Poolevlll996h . Given the probably 
similar masses of warm and cold phase (see Section 12.3)) , it 
appears possible that the latter cooled and condensed from 
the former. This seems difficult for sources in nearby galaxy 
clusters, but may work at higher redshift, where the clus- 
ter gas is expected to be cooler, and therefore to have a 
shorter cooling time. Weak shocks caused by the jet may 
further lower the cooling time. This idea would explain the 
dependence of the appearance of cold gas on redshift (see 
Section 1231). 

A better understanding of these inter-relations would 
not only advance radio source physics, but bear the poten- 
tial of using jets and the associated emission as probes of 
the surrounding (cluster) gas. As cluster gas properties at 
redshifts of one and beyond are s till hard to determine with 
present day X-ray telescopes fe.g. lAndreon et al.ll2005r h this 
may give important constraints, with implications for cos- 
mology and galaxy formation. 

We plan to address the problem in three ways. Here 
we concentrate on simulations that contain the three phases 
from the beginning and study their interaction in kpc-sized 
box simulations. Clues from large scale jet propagation and 
investigation of the possibility of condensing the cold phase 
entirely out of the warm one are deferred to future publica- 
tions. 

The three phases are reviewed individually in Section [2] 
Section [3] discusses the relation to turbulence research, the 
simulations are presented in Section 2] and discussed in Sec- 
tion [5] 



2 THE THREE PHASES PRESENT IN RADIO 
GALAXIES 

2.1 Hot magnetised radio plasma 

Many of the more powerfu l FR II sources ([Fanaroff &; Rilevl 
Il974l : iGilbert et all l2004h possess fat radio cocoons. Pro- 
vided the jets are light compared to the ambient gas, beam 
plasma that has passed the termin al shock region assembles 
in roughly cylindrical cocoons fe.g. lNorman et al . 1982). Ac- 
cording to momentum balance, the head advance speed 
drops with the beam density />j as: 

Vh Vrjevj , (1) 

where 77 = Pj/po, po is the ambient density, v$ is the 
jet speed, and e the model dependent thrust distribution 
factor (< 1). Therefore, the lighter the jets, the slower 
the head advance speed, and the wider the cocoons have 
to be in order to store the shocked beam plasma. Sim- 
ulations still struggle to reproduce the observed cocoon 
widths, which at least in some cases require 77 < 10 ~ 4 
(Krause 2005) , in agreement wi th analytical considerations 
(e.g. Alexand er &; Poolevl fl996). The jet velocity must be 
a reaso nable fraction of the speed of light - jet promi- 
nence (Hardcastl e et al.l 1 19991) and arm- length asymme- 
tries ( Arshaki an &; Lon gair 2004) suggest Vj > c/2 whereas 
a bulk Lorentz facto r of about 10 is su ggested by X- 
ray-Compton studies (|Tavecchio et al . 2004, and references 
therein) a nd the morph ology and size of the bow shock in 
Cygnus A (|Krausell2005h . These conditions in the jet imply, 
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Figure 1. Employed cooling curve. The low temperature part is 
density dependent. 

after the jet shock, a very hot tenuous cocoon plasma with 
an effective temperature of order 

Tc = 3 x 10 11 f^j-) K . (2) 



2.2 Warm X-ray gas 

The radio cocoons are in rough pressure equilibrium with 
the surrounding shocked X-ray gas at, typically, 10 7 -10 8 K. 
Therefore, and for dynamical reasons (compare above) 
the typical density ratio between cocoon plasma and sur- 
rounding X-ray gas should be about 10 -4 . The dividing 
contact surface is subject to Ray leigh- Taylor and Kelvin- 
Helmholtz in stabilities, although it may b e stabilised by 
decele ration j Alexande r! 12002k iKrausd I2005T ) and magnetic 
fields ((Krause &; Camenzindl 120011 ) for some time. Hydro- 
dynamic simulation indicates that a few percent of the 
gas mass pushed aside by the radio cocoon is entrained 
back into the cocoon region by these instabilities, typ- 
ically of the order 10 9 Mq for sources of size 100 kpc 
(Krause 2005) . This picture is confirmed by radio obser- 
vations fe.g. IGilbert et alJl2004h . While the regions close to 
the hotspots of FR II sources are usually sharply bound, 
indicating a stable contact surface, towards the centres the 
emission diffuses away and gets filamentary, indicating en- 
trainment of the surrounding shocked ambient gas. The lat- 
ter is very obvious f rom the high dyna mic range observa- 
tions of Cygnus A of IPerlev et al.1 (Il984h . The entrained X- 
ray gas has been resolved by Chandra observations of the 
same source ([Smith et alJ l2002). In the following, this will be 
called the warm phase in contrast to the hot, radio-emitting 
cocoon-plasma. 

2.3 Cold emission line gas 

A third gas com ponent is reveal ed by the alignment effect 
(for a review see McCarthy 1993). Beyond a redshift of 0.6, 
optical emission is aligned with the radio jet and often cospa- 
tial with the cocoon. This emission is indicative of compara- 
tively cold (T w 10 4 K) gas. The emission is partly polarised 
to varying degree, which is interpreted as a scattering con- 
tribution of light from a hidden quasar by dust grains (e.g. 
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ISolorzano-Inarrea et al.ll2004l, and references th erein). Addi- 
tionally, nebular emission (llnskip et al.ll2002allbh and newly 
formed stars (e.g. lBest et al.lll996l : Dev et alJll997f) may con- 
tribu te to the continuum light (but see also iNeeser et all 
1999). The emission line luminosity is well correlated with 
the radio power. Such a description sounds much like the 
ordinary inter-stellar medium (ISM) in galaxies. The ISM 
is traditionally regarded as a multi-phase medium itself. 
However, recent advances in turbulence research have re- 
vealed rapid transf ormation between cl o uds a nd inter-cloud 
medium, wherefore Elmegre en &; Scald (2004) drop the dis- 
tinction. We adopt this viewpoint and regard all the gas that 
is able to cool significantly by emission on a reasonable time- 
scale (see below) as the third phase, which will be referred 
to a s the cold phase. 

iBaum fe McCarthy! (|2000h give a range of 1O 8 -1O 1O M 
for the cold gas component for radio galaxies at z « 1. This 
is likely to be much more than the radio plasma by a factor of 
at least a hundred, given the density ratio discussed above, 
and about the same as the entrained X-ray gas, as inferred 
from simulation. 



3 CONNECTION TO TURBULENCE THEORY 

'The appearance of instability in viscous shear flows usu- 
ally re sults in th e nonlinear outcome: the onset of turbu- 
lence' (|Shulll992L p 119). FR I jets have been recognised 
as tu rbulent flows early on (e.g. iBicknelll 1 1 984l ; iKomissarovl 
1990). The cocoons of FR II jets are even more likely to be 
turbulent, since these shear flows probably have low veloc- 
ities with respect to the surrounding shocked ambient gas. 
This has already been mentione d in the early simulation pa- 
pers (e.g. IKossI &; Muilerl ll988). The appearance of smaller 
and smaller vortices with increasing resolution has actually 
been o bserved in a resolution study bv lKrause &; Camenzindl 
( 2001), who varied the resolution of their simulations by up 
to a factor of 40. Hence, the turbulent nature of the cocoons 
of FR II sources is evid ent, and the recent literature treats it 
as a m atter of course (|Saxton et al. 2002; Carvalho &; O'Deal 
l2002h . 

An e xcellent review of i nterste llar turbulence can be 
found in [Elmegreen Scalol (|2004f ). They point out that 
compressible turbulence, which is relevant here, is still not 
fully understood. For example, an inertial range where ki- 
netic energy is transferred in a scale-free manner (the Kol- 
mogorov law) to smaller scales (or larger scales as may be the 
case for 2D turbulence) is not expected a priori. The char- 
acter of turbulence is furthermore dependent on magnetic 
fields and self-gravity. Since turbulence is a phenomenon 
which covers a wide range of physical scales, resolution in 
numerical simulations is an issue. In this paper we present 
simulations which while, by necessity, are of limited reso- 
lution, we believe offer improved physical insight. Some as- 
pects of the p roblem we address have been studied by other 
authors (e.g. Kri tsuk fe Norma n 2002, 2004; Mellema et al. 
2002), thereby allowing some consistency check on our re- 
sults. 

Turbulence quickly decays, if not driven constantly. Its 
properties, however, may depend on the driving mechanism 
( Elm egreen &; Scaloll2004f ) . Many turbulence simulations use 
large scale driving by source terms. In the simulations pre- 
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Figure 2. Sketch of the setup for all simulations. The three 
phases soon mix due to fluid instabilities, and form a turbulent 
region powered by the initial gravitational and kinetic energy. 



sented below, the energy is initially stored mainly in ki- 
netic energy, and transferred slowly to turbulent motions by 
the Kelvin-Helmholtz instability. Because the linear growth 
times of these instabilities depend on the length scale (Equa- 
tions J6}, the driving starts on small scales and involves 
larger scales the longer the simulation is run for. This is de- 
signed to reproduce the conditions in jet cocoons as closely 
as possibly, but it is different to most simulations in the 
literature. 

There are a number of general results concerning com- 
pressible turbulence from theory and simulation which we 
expect to find reflected in our simulations. Turbulence the- 
ory finds log-normal density probabilit y distribution func- 
tions (PDF) for isothermal turbulence ([Elmegreen &; Scalol 
2004): power law tails are expected for specific heat ratios 
other than one. The volume filling factor as a function of 
density is determined to be a power law from simulation. 
The velocity PDF is expected to be Gaussian from analytic 
work but is found to be more exponential from simulations. 
In the cold phase, t he average Mach number grows with 
density as M oc p 1/2 (Kri tsuk fe NormadlioO? ). 



4 THE SIMULATIONS 

Hydrodynamic simulations were performed in order to study 
the turbulent interaction of the three phases. The equations 
for mass (p dx 3 ), momentum (pv dx 3 ), and internal energy 



(e dx ) are 
t + V.( P v, 



<9pv 



+ V(pvv) — — Vp — gp 



de 

m 



+ V(ev) = -pV • v - 



(3) 
(4) 
(5) 
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where m p is the proton mass and the system is closed by the 
equation of state p = (7 — l)e with 7 = 5/3 being the ratio 
of specific heats. These equations were evol ved in two Carte- 
sian dimensions with the code NIRVANA ([Ziegler &; YorkeJ 
Il997h . This code is second order accurate in space and time, 
and is able to resolve shocks via artificial viscosity (com- 
pressible hydrodynamics) . 

Cooling is treated as a source term in the energy 
equation. The cooling curve was adopted from Basson 
(l2002h. For temperatures b etween 10 4 K and 10 8 5 K, the 
ISutherland &; Dopital (|l993h cooling curve for solar metal- 
licity was used. Above that temperature the cooling is 
due to bremsstrahlung: A = 2.05 x 10 -40 a/T(1 + 4.4 x 
10 -10 T) Jm 3 /s, where the second term in the bracket is 
the relativist ic correction. For one run, the cooling curve 
was cut below 10 4 K, to emulate effects of photoionisation. 
For the others, cooling by collis i onal d e-excitation of molec- 
ular hydrogen (jTegmark et al.1 Il997t 10 2 K-10 4 K) and by 
emiss ion lines from H2, HD, and CO molecules (|Puv et al.l 
1999, < 100K) were included. This cooling curve is shown 
in Figure [1] Another source term is the constant gravity of 
g = 10 -10 ms -2 . The gravity does not influence the simula- 
tions by much but was included to make the simulations as 
realistic as possible. 

We have also used the third order accurate Flash code, 
which is forced to conserve energy and momentum to check 
the results. We encounter numerical problems with both 
codes, when regions of hugely differing densities approach 
each other. This restricts the density ratios in the simu- 
lations. It turned out that, with given computational re- 
sources, we could perform better simulations with Nirvana. 
We could not reach enough resolution with Flash to produce 
reliable results in the parameter range of the presented Nir- 
vana simulations. Therefore, we discuss mainly our Nirvana 
simulations. 

4.1 Setup 

The initial setup of the simulations is shown in Figure [2] 
The idea is to model the situation at the contact surface be- 
tween jet cocoon and shocked ambient gas. The upper half 
of the computational square box of size [1 kpc x 1 kpc] , re- 
solved by [512x512] cells, is initially filled with gas of density 
10 4 ra p m -3 . The pressure was set to decrease linearly with 
height, in order to achieve hydrostatic equilibrium. With a 
temperature of 2 x 10 6 K, the cooling time is about 10 Myr, 
which is similar to our simulation time. Five elliptical clouds 
with a very high density (see Table [TJ) and various radii 
(0.008,0.02,0.04,0.07,0.004 kpc) have been placed in the high 
density region at a height of 0.6 kpc. This results in a total 
cold mass of 10 4 -10 7 Mq kpc -3 . For determination of this 
density and similar ones to follow, the cells are assumed to 
represent cubes of three equal sides. Hence, the values are 
upper limits due to the unknown 3D structure. Observed 
emission line gas masses (see above) and radio cocoon sizes 
of about Vcoc = Kr 2 coc h « tt(15 kpc) 2 100 kpc ^ 3 x 10 5 kpc 3 
lead to cold gas masses of 10 3 -10 5 Mq kpc -3 . We also in- 
clude a control run without any clouds. 

The lower half of the box contained low density gas 
(m p m -3 ) at about the same pressure, which was adjusted 
experimentally for best numerical stability of the initial con- 
figuration. This gas was given an initial velocity of v = 



2 x 10 7 ms -1 , corresponding to a Mach number of 0.8 (80) 
with respect to the gas in the lower (upper) part of the 
grid. The boundary is perturbed as a sine wave of ampli- 
tude 0.03 kpc. This excites Kelvin-Helmholtz instabilities, 
with a linear growth time of: 

^= 5Myr (l^) (iF^ GxJms-O " 1 ' < 6 > 

where 77 is the density ratio, and A is the wavelength of the 
perturbation. The perturbations are seeded with a wave- 
length of the box size. However the discretisation seeds per- 
turbations on the scale of the resolution, which dominates 
the evolution. The initial velocity is in the horizontal direc- 
tion only, not parallel to the surface. As the smallest pertur- 
bations grow fastest, this speeds up the development towards 
turbulence. The gravitational acceleration is comparatively 
low, wherefore the linear Rayleigh-Taylor instability does 
not contribute significantly to the driving. 

Simulation parameters that were varied between simu- 
lations are given in Table [1] Note that for the model with 
the temperature cut in the cooling function the cloud den- 
sity was set to yield the minimum temperature. 

The timestep was found to decrease substantially 
with contained mass. Hence the high mass simulation was 
stopped after 3 Myr. 

4.2 General dynamics 

We begin by discussing run 7c4 which serves to illustrate 
the dynamics shown in all simulations containing clouds (see 
Figure[3]). At the beginning of the simulations, the dynamics 
is governed by the Kelvin-Helmholtz instability. There is also 
resonant a kink instability which could contribute in prin- 
ciple. This instability produces a ch aracteristic shock pat- 
tern (a kinked shock) at the interface ([Bassett &; Woodward! 
1995). However, we have checked images of the velocity di- 
vergence, and could find no trace of this pattern. The evolv- 
ing Kelvin-Helmholtz instability produces vortices in the 
subsonic lower half of the computational domain, and sends 
strong shocks into the upper medium and the clouds. Si- 
multaneously, the clouds cool, either to the temperature cut 
of 10 4 K, or to about 10 3 K. In this phase, the clouds evolve 
m uch like the ones in t he shocked cooling cloud simulations 
of lMellema et al. l (120021). Once t he sho ck passes, the clouds 
are compressed. iMellema et al.l (|2002n then see a fragmen- 
tation into small stable cloudlets: our simulations do not 
have sufficient resolution to observe fragmentation. Instead 
we observe the formation of one filament per cloud. When 
the contact surface reaches it, the filament spreads along 
the surface. The filament is initially unable to penetrate the 
contact surface due to the long Rayleigh-Taylor time-scale 
(Figure [3] second row, left). Subsequently, one of the clouds 
is stretched and pushed towards the upper boundary by a 
rising part of the instability. This particular feature is espe- 
cially related to the numerical issues detailed in section l4.10l 
Another cloud, located at a declining part of the contact 
surface is reached by a branch of the low density gas and 
drawn into it (Figure[3] second row, right). There it is greatly 
stretched and soon forms a filamentary system. When the 
first of the dense filaments reaches the lower boundary (be- 
tween the centre and the right plots of the bottom row in 
Figure [3} the structure of the flow changes. The low density 
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Figure 3. Logarithmic density plots at nine representative simulation times (0, 0.2, 0.9, 1.7, 2.4, 3.2, 4.1, 4.8, and 7.7 Myr) for RUN 7c4. 



Table 1. Simulation parameters 



Label 


Cloud 


Total 


Upper 


Temperature 


Final 




density 


cloud mass 


boundary 


cut 


time 




[ra p m -3 ] 


[M kpc- 3 ] 




[K] 


[Myr] 


ctrP 


10 4 


1.5 x 10 3 


closed 





10 


5™ 


10 5 


1.5 x 10 4 


closed 





10 


6o0 


10 6 


1.5 x 10 5 


open 





10 


7c4 


5 x 10 6 


7.4 x 10 5 


closed 


10 4 


10 


7o0 


10 7 


1.5 x 10 6 


open 





10 


7c0 


10 7 


1.5 x 10 6 


closed 





6.4 


8o0 


10 8 


1.5 x 10 7 


open 





3 



a contains no clouds, just background density 

gas can no longer rush through from left to right, and the 4.3 Dynamics of the different setups 
motions turn completely into turbulence. Towards the end, 

the upper part moves left and the lower one to the right, Logarithmic density plots at 3 Myr are shown in Figure H 
forming one large scale vortex. f° r an six simulations containing clouds. The mass loading 

clearly affects the evolution of the Kelvin-Helmholtz insta- 
bility. This is most evident when comparing the right column 
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Figure 4. Logarithmic density plots at 3 Myr for all the simula- 
tions with clouds in the order of Table [T] (5c0: top left, 6o0: top 
right, 7c4: middle left, 7o0: middle right, 7c0 bottom left, 8o0: 
bottom right). An open upper boundary (right plots) amplifies 
the Kelvin-Helmholtz instability. More cold mass (towards the 
bottom) speeds up the formation of filaments. 



of Figure U (open upper boundary). When the first filament 
is stretched towards the upper boundary, the low density gas 
flows along, leaving the grid if the upper boundary is open. 
In the simulations with lower cloud density, this feature is 
suppressed, and the low density gas stays in the lower part 
of the computational domain. Also for the simulations with 
lower initial cloud mass the overall evolution is slower. 

The direct comparison between open and closed bound- 
ary (runs 7c0 and 7o0) shows that dense gas assembles near 
the boundary for the closed simulations and that this gas 
leaves the grid for the open ones. This is not compensated 
by inflows, and hence, the open boundary cases have less 
cold mass than the closed upper boundary simulations. 




Figure 5. Two-dimensional velocity power spectra for run 7c4 
at 0.2 Myr (top left), 3.2 Myr (top right), and 7.7 Myr (bottom 
left) and for run 7o0 at t = 7.7 Myr (bottom right). The power 
spectra are evaluated in the zero momentum frame. Wave-vector 
zero is at the centre. 
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Figure 6. Energy spectrum (E(k r )dk r = 
J P(k r ,k ( f ) )k r dk ( f ) dk r ), X- and Y-cuts through the power 

spectrum of run 7c4 at 7.7 Myr. The large scales isotropise first, 
showing a Kolmogorov-like decline (k~ 8 / 3 in the power spectrum 
cuts, corresponding to fc _5//3 in the energy spectrum). The 
smallest scales do not isotropise up to the end of the simulation, 
which is due to the bar feature. 



The temperature cut results in a significantly reduced 
maximum density (by a factor of order 10). The dens- 
est cloud occupies a bigger area, i.e. it is less condensed. 
Also, the aforementioned uplifting of low density gas is sup- 
pressed. 
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Figure 7. Root mean squared Mach number versus time for all 
runs with clouds. The Mach number generally grows with time. 
The simulations with higher mass load have higher average Mach 
numbers. The high amplitude features are due high density fil- 
aments hitting a boundary, and are especially pronounced for 
closed boundaries.. 



4.4 Power spectrum 

For comparison to turbulence studies, it is important to 
check isotropy and scale-dependence. Therefore, we exam- 
ined the power spectrum of the velocity, 



p(k) 



• „. „, J bo 



Vi e lkx d 2 x\ 2 . 



We use velocities in the zero momentum frame although 
the results do not depend on that particular choice. Two- 
dimensional representations of the power spectrum are 
shown in Figure [5] The behaviour for the different runs is 
very similar. Hence, only results from runs 7c4, and 7o0 are 
shown. Even the early power spectra look quite isotropic, 
apart from a horizontal bar feature. However, the energy 
spectra up to a few million years are dominated by the bar 
feature, which produces power law indices around —2. Large 
scales isotropise first, and at 7.7 Myr, we find a Kolmogorov- 
like isotropic power spectrum for most scales (Figure [6}. 
The bar feature varies in strength, but is generally decreas- 
ing. The open upper boundary simulation showed a stronger 
feature at the same time. The same is true for simulations 
with less cold mass load. 

We conclude that for larger scales, isotropic, 
Kolmogorov-like turbulence is essentially reached after 
a few million years of simulation time. Some anisotropy 
remains on small scales. 



4.5 Mach number, density and pressure 
distributions 

The root mean squared (rms) Mach number is shown in 
Figure [7] It is generally well above one with the exception 
of run 5c0. 

The choice of boundary condition has a big impact on 
the rms Mach number. We compare the three simulations 
with roughly the same initial density (7o0,7c0,7c4). The big 
bump in the 7c0 line is produced by one large filament hit- 
ting the upper boundary. The latter releases so much kinetic 
energy that the rms Mach number reaches up to a factor of 



ten larger than in the open boundary case for the same den- 
sity. Remarkably, after relaxation, the two simulations rejoin 
the same path. The bump is completely damped in the sim- 
ulation with the temperature cut (7c4). The temperature 
cut reduces the peak values of the density, which leads to 
the more modest interaction with the grid boundary. 

The most evident feature in Figure is the constant 
rising of the rms Mach number and the dependence on the 
mass loading. More cold mass results in higher rms Mach 
numbers. This observation can be explained by a remarkable 
relationship between the density and the Mach number. 

The relationship between density and Mach number is 
shown in the 2D-histograms of Figure [8] These diagrams 
display the volume occupation of gas at a given density and 
Mach number. We have chosen the time of 3 Myr (10 Myr) 
for the higher (lower) cold mass load simulations, because 
the turbulence is well evolved at that time, with the mix- 
ing between the low and medium density gas still being ac- 
ceptably small. These plots look very similar at most other 
simulation times. Mixing of these gas phases is indicated 
by the two peaks that are initially at 10 -27 kgm -3 and 
10 -23 kgm -3 coming closer together. 

Between these densities, the Mach number is indepen- 
dent of the density. Since the simulation is effectively iso- 
baric at those densities (see Figure [9} , this means that the 
kinetic energy density is roughly constant, indicative of a 
quasi-equilibrium state. The same relation would be ex- 
pected in a shock dominated scenario, where T oc v 2 . At 
high densities, the internal energy is reduced by the radi- 
ation losses, which increases the Mach number. Although 
this can change the pressure by many orders of magnitude 
in individual cells, the bulk of the material can compensate 
for the pressure loss by contraction and shock heating, which 
brings the gas back to near isobaric conditions. Since the ve- 
locity is not correlated with the density, the dependence of 
the Mach number on density (M 2 = pv 2 /*yp) is dominated 
by the explicit density depe ndence. These findi ngs agree well 
with the 3D simulations of iKritsuk fc Normanl ([§004) who 
drive the turbulence by the thermal instability, only, and do 
not include the low density phase. 

The M oc p 1 / 2 branch has got a higher occupation and 
extends to higher Mach number for higher mass loading 
and later simulation times. This causes the behaviour of the 
rms Mach number in Figure [7] 

Typical pressure versus density histograms are shown 
in Figure [9] Although the occupied phase space spans many 
orders of magnitude in both dimensions, the pressure is re- 
markably peaked towards a central value, even in regions 
with high density and a shorter cooling time - here the pres- 
sure is, in general, at most a factor of ten below the low den- 
sity gas. The prominent linear features in these figures are 
due to two processes. One is quasi-adiabatic expansion and 
compression of the marginally cooling gas at 10 -23 kgm -3 , 
which initially fills the most of the upper part of the grid. 
The other is the quasi-isothermal regime at the prominent 
drop of the cooling function at about 10 4 K. The line is 
most prominent in the simulation with the cut in the cool- 
ing function but is also apparent in the other simulations at 
a temperature of about 14,000 K. Evidently, at this point the 
energy gain due to shocks caused by the turbulent motions 
balance the radiation losses. 

The density PDFs are dominated by a power law be- 
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Figure 8. Phase diagrams showing the logarithmic volume occupation at a given density and Mach number for runs 5c0 and 6o0 (top 
at 10 Myr), and runs 7c4 and 7c0 (bottom at 3 Myr). The distribution is bimodal: low density gas is unable to cool and forms a 
turbulent system where the Mach number distribution is independent of the density. Gas with short cooling times joins at a density of 
about 10 -23 kgm -3 and its Mach number grows roughly with the square root of the density. The fiducial line indicates a M oc p 1 / 2 
dependence. Colour coding is the same as on other figures. 



haviour (Figure [TO]) . As long as the peaks of the initial con- 
dition remain visible, the region between the peaks tends to 
follow a ]/(/)) oc p -1 ^ 2 law, where V{p) is the occupied vol- 
ume at a particular density p. The high density part of the 
distribution reflects the mass loading, and is independent 
of the low density part. This region also follows the p _1//2 
law, at least for higher mass load. Over time, the region 
around 10 -22 kgm -3 becomes more and more populated. 
As the peaks diffuse, the distribution becomes uniform with 
exponential type cutoffs towards low and high density. 



4.6 Temperature distribution 

We show the distribution of gas mass as a function of tem- 
perature in Figure 1111 The general behaviour is well de- 
scribed by a dM/dT oc 1/VT law. The high temperature 
part (T > 10 6 K) evolves independently from the cooler 
part. Due to numerical mixing, the upper temperature of the 
gas moves to lower temperature keeping an exponential-type 
cutoff. In all simulations, the distribution evolves quickly to 
follow a dM/dT oc 1/VT form for T > 10 6 K up to the 
cutoff. There is a well defined peak at T « 10 6 K which 
broadens during the simulation. For T < 10 6 K the distribu- 
tion evolves more slowly towards the same dM/dT oc 1/VT 
established at higher temperatures and at a rate which in- 



creases with initial cloud mass. A peak at 14,000 K is well- 
defined in all simulations, and does not evolve with time. 
Below ~ 10 3 K the distribution evolves to form what may 
be a power- law decline at low temperatures. 



4.7 Filling factor and mass 

NIRVANA conserves the total mass to machine accuracy. In 
the simulations with closed boundaries, the mass is therefore 
constant (Figure IT2l right). The other simulations lose up 
to about 50 per cent of the mass over the grid boundary - 
the ordering according to initial mass loading is essentially 
maintained throughout the simulation. 

Since radiation transfer is not included into the simula- 
tion, the ionisation structure and nebular emission proper- 
ties cannot be deduced in detail. However, from the strong 
drop of the gas mass towards higher temperature (Fig- 
ure [TT]), and the drop of the cooling function towards lower 
temperature, it is clear that the gas with a temperature close 
to 14,000 K is responsible for much of the optical luminos- 
ity. In Figure [12] (left) we show the volume filling factor of 
gas with temperature 14,000±2,000 K, which corresponds 
to the peak in the mass versus temperature distribution in 
Figure [TT] This filling factor depends mostly on mass, with 
the three most accurate simulations (5c0, 6o0, and 7c4, see 
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Figure 9. Phase diagrams showing the logarithmic volume occupation at a given density and pressure for runs 5c0 and 6o0 (top at 
10 Myr), and runs 7c4 and 7c0 (bottom at 3 Myr). The pressure of much of the simulated volume is within a factor of a few, even in 
the strongly radiating part. Fiducial lines of slope 5/3 and 1 highlight adiabatic expansion features, prominently around a density of 
10 -23 kgm -3 (upper background gas), and isothermal components, respectively. In run 7c4, the 10 4 K line marks the border of the 
accessible phase space, in the other simulations, gas assembles at a slightly higher temperature due to the drop of the cooling function. 
Colour coding is the same as on other figures. 



Section 14.101 below) following nearly parallel paths. As ex- 
pected, higher initial mass loading results in a higher filling 
factor. For most of the simulation times, the filling factor 
stays between « 10 -3 and ~ 10 -2 . Since these are 2D sim- 
ulations, this would correspond to a volume filling factor of 
between 10~ 5 and 10~ 3 in 3D, assuming similar structuring 
in all dimensions. 



4.8 Energy loss 

Energy conservation is not guaranteed formally by the nu- 
merical method employed by NIRVANA. However, for prac- 
tical considerations, it is usually conserved to reasonable 
accuracy, given sufficient resolution. This is the case for at 
least three of our runs (ctrl, 5c0, and 6o0, see Section 14.101 
below for further discussion) - we therefore first focus on the 
energy evolution of these three simulations. Figure [13] (top 
left) shows the time evolution of the total energy for all runs, 
including the control run. Since there are no energy sources 
in the simulated volume, the total energy should monoton- 
ically decrease, at least in the closed box cases, due to the 
radiation term. In the control run without clouds the sys- 
tem loses energy at a nearly constant rate of 4 x 10 28 W. For 
the simulation with the lightest clouds, the rate increases to 



18 x 10 28 W, and for run 6o0, the total cooling increases 
further to 22 x 10 28 W. As might have been expected, the 
presence of the cool clouds increases the cooling rate. This 
is no longer evident for the simulations with higher mass 
loading. While at later times they turn to monotonic energy 
decrease at an even greater rate, they are obviously domi- 
nated by numerical effects for some and in the case of run 8o0 
even much of the simulation. These numerical errors can be 
traced back to the kinetic energy. While the thermal energy 
(Figure [T3] bottom left) shows essentially a smooth energy 
decline, with the only significant dependence being on mass 
in the expected way, the kinetic energy displays strong and 
artificial wiggles, which is discussed below. The extraction 
of energy from the high temperature system can be seen in 
yet another way. The fraction of energy in gas hotter than 
10 6 K is shown in Figure [13] (top right). The system not only 
loses energy faster with increasing mass loading, but addi- 
tionally energy is increasingly found in the low temperature 
gas. 



4.9 Cold mass dropout 

The increased system cooling rate due to the cold mass con- 
tent leads to amplified mass dropout - this is illustrated in 
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Figure 10. Left: Density probability distribution functions at a representative time of 2.5 Myr for all simulations. The fiducial line is 
proportional to p -1 / 2 . Right: The same for run 7c4 for a series of simulation times. 
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Figure 12. Filling factor (left) of gas in the 14,000 K peak and total mass (right) over time for all simulations with clouds. The mass is 
conserved by NIRVANA, and therefore, only the simulations with open boundary lose some mass (up to « 50%) 
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Figure 13. Energy evolution over time for all runs. Top left: total (i.e. thermal, kinetic and gravitational) energy; top right: fraction of 
total energy in gas hotter than 10 6 K; bottom left: thermal energy; bottom right: kinetic energy. All runs besides ctrl, 5c0 and 6o0 show 
problems with energy conservation (see text), which can be traced to the kinetic energy. 
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Figure 14. Time evolution of the cold gas phase (T < 10 6 K) for cloud containing simulations. Left: total amount of cold gas for the 
simulations with closed boundary. Right: Evolution of the cold gas mass fraction for all simulations. In the control run, no cold gas was 
detected over 10 Myr. 
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Figure fl4l As discussed in Section 14.71 the total mass is con- 
served in the closed boundary simulations. The cold mass is 
therefore increased by dropout from the warm phase, only. 
In the control run, there is no cold mass present initially, and 
none drops out during the simulation time. Runs 5c0, 7c4, 
and 7c0 start with a cold mass of 30, 1434, and 2957 Mq, 
respectively: they gain cold mass at rates of 13.4, 23.3, and 
27.3 Mq Myr -1 , on average. The higher mass loaded simu- 
lations are limited by the available mass. At the beginning 
of the simulation, their cold gas mass is already ~ 90 per 
cent of the total mass. They gain almost all of the rest dur- 
ing the simulation time. The cold mass evolution for run 5c0 
can be well fitted by an exponential growth with e-folding 
time of 6 Myr. This corresponds to a simple interpretation in 
which the rate of mass dropout is proportional to the mass 
of cold gas. We note that runs 7c0 and 7o0, which differ 
only in their boundary condition, have very similar cold gas 
mass fraction histories, suggesting that the boundary condi- 
tions are unimportant for the fractional mass dropout. This 
is confirmed by the similar shape of all the curves in Fig- 
ure [14] (right). The simulations with small initial cold gas 
mass fraction show increasing growth up to about 80 per 
cent. Further growth proceeds in a linear fashion. 

4.10 Numerical issues 

As indicated above in Section 14.81 the higher mass loaded 
simulations show increasing problems with energy conserva- 
tion, in particular kinetic energy. Careful inspection of the 
simulations showed that the velocity errors depend strongly 
on the density contrast between cells. In the high mass 
loaded simulations, as soon as the very light gas in the lower 
half of the simulations comes into contact with the high den- 
sity clouds, high velocity errors result. This is most evident 
in the total and kinetic energy plots (Figure [13]). The higher 
the initial density, the higher the numerical error. Simula- 
tions 5c0 and 6o0 show no indication of the problem. Run 
7c4 shows moderate errors. Runs 7o0 and 7c0 display moder- 
ate errors with a singular spurious event around 2 Myr. The 
fact that 7c4 shows only moderate errors at this time deter- 
mines the critical density contrast to be about 10 7 . Run 8o0 
is dominated strongly by numerical errors for much of the 
simulation time. We have checked that these errors decrease, 
as expected, if the numerical resolution is increased. Some 
aspects of the simulations are affected worse than others. Re- 
garding most of the results discussed so far, a noticeable sim- 
ilarity and continuation of the results argues for the kinetic 
energy errors not to have a too large effect. In particular, 
the power spectrum is completely unaffected. Comparison 
of the phase diagrams for runs 7o0 and 7c0 shows that once 
the gas which was accelerated in the singular spurious event 
has left the grid (7o0), the diagram soon reverts to normal, 
in contrast to the closed boundary simulation (7c0). The 
high mass load simulations give the best signal to noise in 
the phase diagrams. However, because of the energy con- 
servation issue, we include run 8o0 only to demonstrate the 
limitation of the method most clearly, and base our conclu- 
sions exclusively on the very best results. 

A further check we have perform ed was to repeat 
the simulation wit h the FLASH code ([Frvxell et al.l boQQ: 
I C alder et al. I l2002h , which guarantees energy conservation. 
However, in this case another problem appeared in low reso- 



lution runs. Here, the clouds were soon surrounded by a low 
pressure shell of a width equal to the resolution limit, which 
seriously damped the propagation of shocks into the clouds. 
The clouds essentially kept their shape for long times, even 
if accelerated. This problem could be solved by going to very 
high resolution, which corresponds for regions of smooth flow 
roughly to a four times better resolution compared to the 
simulations presented above, due to the higher order scheme 
of FLASH. We also adjusted the setup to have a single more 
massive cloud and allowed for higher background density 
(top: 10 5 m p m -3 , bottom: 10 2 m p m -3 ). With such a setup 
we obtain very similar results to the NIRVANA simulations 
(see Figure [T5]) . The cloud soon collapses to form a filament, 
dispersing into little cloudlets afterwards. The Mach-number 
density relation is consistent with the NIRVANA results, 
and the gas mass versus temperature distribution shows the 
same general behaviour and structure. 

4.11 Velocity distribution functions 

The velocity PDFs for the two most reliable runs (5c0 and 
6o0) are shown in Figure 1161 As expected from general tur- 
bulence theory, they show a Gaussian core with signs of 
acceleration due to drag and upwards driven shock waves 
towards positive velocities. These two most trustworthy sim- 
ulations also show a dependence of the width on the mass 
loading. This confirms the earlier findings that in the simu- 
lations with higher mass load, the fraction of kinetic energy 
in the cold phase rises faster. 



5 DISCUSSION 

We present simulations of the Kelvin-Helmhotz instability 
with clouds of differing density. The clouds are shocked, col- 
lapse into a filament, and then disperse into cloudlets and 
more filaments. There is no indication of the cold mass be- 
ing heated beyond typical emission line gas temperatures. 
On the contrary, over time more and more gas condenses 
into the cold phase. 

Scale-free turbulence is establishes on the larger scales 
first. The smallest scales remain anisotropic up to the end 
of the simulation 

In general, the rms Mach number increases during the 
simulation. This is due to the cold gas aquiring high veloci- 
ties. In particular, for gas with longer cooling time, the Mach 
number is independent of the density, whereas for strongly 
cooling gas, the Mach number is proportional to the square 
root of the de nsity. The latter dependen cy is identical to the 
one found bv lKritsuk & Norman] (|2004h . 

The density probability distribution function shows a 
power law behaviour between the peaks imposed by the 
initial distribution, P P (p) oc p -1 ^ 2 . A similar behaviour is 
seen in the gas mass versus temperature plots. Generally, 
the distributions follow M(T) oc T _1//2 , with the region be- 
tween 10 4 K and 10 6 K evolving gradually over a few million 
years to follow the same scaling as at high temperature. But 
what is the cause of this behaviour? Because of the roughly 
constant pressure, the entropy distribution is also a power 
law. This rules out adiabatic processes. We are left with 
shock heating, cooling and mixing. Cooling and shocks con- 
tribute little for the gas with T > 10 6 K. So, mixing seems 
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Figure 15. Selection of results for the validation run with 
FLASH. Top: logarithmic density distribution at 1 Myr, middle: 
density Mach number histogram at 1 Myr, bottom: gas mass over 
temperature, all gas, bottom right: gas mass over temperature 
between 5,000 K and 30,000 K, staggered for 13 timesteps from 
0.1 Myr to 1.3 Myr. 



to be the best explanation for the power law slope of the 
distribution functions of density and temperature. Mixing 
is inevitable in multi-phase turbulence simulations. In real- 
ity, mixing is likely to occur at a lower level, consequently 
the enhancements over the mixing distributions should be 
more pronounced. There is a worry that the peak around 
14,000 K described in the following might be entirely due to 
the mixing. Against this possibility argues: 

• The peak builds up early, when the mixing in the cold 
gas is still at low level. 

• The 14,000 K peak joins the lower temperature gas 
above the mixing power law. 

• The much less diffusive Flash simulation, produces a 
very similar peak. 

We conclude that the 14,000 K peak is probably unaffected 
by the mixing. 

This peak is in fact likely to be due to shock heating 
combined with the effects of the steep cooling curve. The 
peak appears in all simulations at most times. Since there 
is a drop in the cooling function towards lower temperature 
and there is a drop in gas mass towards higher temperature, 
the gas close to this temperature may contribute most of the 
optical emission. Its filling factor depends on the initial mass 
loading, and corresponds to 10 -5 to 10 -3 , for a realistic 3D 
generalisation. The velocity of this gas is distributed like a 
Gaussian in run 5o0. Run 6c0 has a pronounced deviation 
on the positive side due to acceleration, and appears unre- 
laxed. The typical width of observed emission lines may be 
understood by the Mach-number density relation. Gas with 
low cooling time was shown to have M oc y/p, i.e. the ve- 
locity is constant. This is valid for gas of a temperature up 
to ~ 10 6 K. For hot gas, the Mach number is constant, i.e 
the velocity drops with p -1 ^ 2 . The hot radio plasma, has 
velocities of typically 10 5 kins" 1 when entering the cocoon, 
which will become increasingly turbulent. If the entrained 
ambient gas is roughly 10,000 times denser, it would be ac- 
celerated to typical velocities of 1000 kins" 1 . If in z ~ 1 
radio galaxies that gas has a temperature of > 10 7 K as for 
nearby sources, then gas with T < 10 6 K would have a typi- 
cal velocity width of a few 100 kms" 1 . However, the typical 
observed velocity width is ~ 1000 kms -1 . Hence, either the 
density ratio in these sources is typically less extreme than 
assumed here, or the typical temperature of the ambient gas 
is below 10 7 K (for both, a factor of ~ 10 would suffice). 

Excess cold mass is located in a broad peak around 
10 3 K, which would contain substantial amounts of molecu- 
lar hydrogen. Most of the gas mass would be at this tem- 
perature, unless there is heating by photoionisation. 

Higher mass loaded simulations extract the energy 
much more efficiently from the warm and hot phase. There- 
fore, the velocity widths as well as the energy extraction 
increase with mass loading. The overall cooling rate of the 
system was increased by about a factor of five for the low 
mass loaded simulations. For the high mass loaded simula- 
tions it was shown to increase further, whereas the quality 
of the simulations gets worse. 

The simulations demonstrate the loss of numerical ac- 
curacy with increasing density contrasts. Three of the sim- 
ulations conserve energy to a very good accuracy; further 
three have moderate problems, whereas the highest mass 
loading simulation has significant errors with respect to the 
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Figure 16. Velocity distribution functions for gas at temperature 14,000 ± 2,000 K for runs 5c0 and 6o0. Left: horizontal velocity. Right: 
vertical velocity. For better signal to noise, output between 1 and 6 Myr (0.1 Myr steps) has been staggered. The shift towards positive 
horizontal velocities, more pronounced for the higher mass load case is the result of the accelaration by the low density gas. 



energy conservation. This is a clear trend. We conclude that 
the highest density contrast in a simulation is the limiting 
factor for numerical reliability. Despite these problems at 
extreme density ratios we believe the main trends and re- 
sults are very reliable - we have checked this by comparing 
results to a high-resolution simulation performed using the 
FLASH code which is forced to be energy conserving. 

The cold gas clouds effectively act as condensation nu- 
clei. Strong growth of the cold gas mass is observed in all 
simulations which is in contrast to the control run in which 
no gas cools below 10 6 K for the whole simulation time. The 
lowest gas mass run, which is furthest from saturation, gains 
a factor of more than five in cold gas mass within 10 Myr. 
In a real radio cocoon the supply of warm gas would be 
ample, and exponential growth would result. Hence, a cold 
gas mass seed of order 100 Mq would be sufficient to con- 
dense the observed ~ 10 9 Mq in the radio galaxies at z ~ 1 
mentioned in the introduction, provided the conditions are 
similar to those considered here and if the typical source age 
is of order 100 Myr. It is worth noting that w 10 9 Mq is also 
the amount of X-ray gas entrained int o the radio cocoon as 
infered from simulation ([Krause 2005). A consistent picture 
emerges, if most of the entrained X-ray gas cools to the cold 
phase via the help of a small cold-gas seed. As this gas cools, 
most of the radiated energy is not in the X-ray band, but in 
the optical. Since in the real radio cocoon, unlike in the sim- 
ulation here, the most energetic part of the system would be 
the radio emitting hot phase, a tight correlation between ra- 
dio and optica l emission would be expected, which is indeed 
observed (e.g. |Mc£arth3|l993). Such a mechanism also of 
course explains the increased visibility of optically emitting 
regions after the passage of the radio cocoon. 

Due to the large temperature ratios, heat conduction 
is a particular concern. However, it may still be negligible 
here, as we argue in the following. Fir st, the evaporation 
lengthscale ([Boehringer &; Fabian|[l989h is at least a factor 
of thousand below our resolution limit. Second, using Spitzer 
conductivity, we can estimate the energy transfer per volume 
due to heat conduction: 

—T~>3 . 5 / T 1 \ 3-5 

. = V.[^VT]^,10-W m -(^) . 



This leads to a heat conduction timescale of: 

= e_ 10- 13 Jm- 3 10 / T \~ 35 

Thc e ~ 10- 23 Wm- 3 (T/10 6 K) 3 - 5 S U0 6 K7 

Since our simulations represent a dynamical equilibrium, we 
should compare this to the cooling timescale in order to find 
out to what extent heat conduction could contribute. The 
cooling time for the hot, warm and cold gas is of order 10 20 s, 
10 15 s and 10 11 s, respectively. So, heat conduction from the 
hot phase would clearly dominate everything else. However, 
in reality the presence of magnetic fields strongly suppresses 
heat conduction in and from the hot phase. If this was not 
the case, radio cocoons could not exist, since the internal 
energy would have been lost on the light crossing timescale. 
Among the cool gas, the heat conduction timescale is too 
long to be of concern. Therefore, we are left with the heat 
flow from the warm to the cold phase. Here, the heat con- 
duction is about 10 times faster than the radiative cooling. 
For the simulation parameters, saturation sets in just above 
10 6 K (?), and does therefore not influence the conclusion. 
However, also in the warm gas, the heat conduction is sever ly 
suppressed due to the presence of magnetic fields, at least 
for the current cosmological age. Measurements from X-ray 
data su ggest a suppression by one to three orders of mag- 
nitud e (jEttori fe Fabianll2000l : iMarkevitch et al.ll2003l : iNathl 
l2003h . The magnetic field in the warm gas around high red- 
shift radio sources is unknown. However, if its effect on heat 
conduction is similar than in nearby cases, it will probably 
be a minor contribution to the total heat transfer. If the heat 
conduction would be less strongly suppressed, the processes 
described in this paper would be further amplified - more 
heat would be radiated by the cold phase and more warm 
gas would condensate on the cool gas. 

The simulations presented here are 2D, only. This is an 
important first step, and it is very likely that the mecha- 
nisms we discuss would also be present in 3D - the Mach- 
number density relation is identical to published 3D results. 
To get more accurate estimates for example for cooling and 
mass dropout amplification or filling factors, 3D simulations 
are required. This is likely to be challenging, given the im- 
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portance of high resolution, but will no doubt be possible in 
the near future. 



6 CONCLUSIONS 

We presented 2D turbulence simulations initiated and driven 
only by the Kelvin-Helmholtz instability. We add dense, 
cooling clouds to the simulations, so that the interaction of 
the three phases present in the cocoons of high redshift radio 
galaxies may be examined. We find a turbulent cascade with 
energy spectrum proportional to k~ 5 ^ 3 at larger scales, as 
in Kolmogorov turbulence. Mach number and density show 
a bimodal correlation with a break at a characteristic den- 
sity, where the cooling time is about the simulation time. 
In the high density part, the Mach number goes as p 1 ^ 2 , 
whic h is the same as found b y independent 3D simulations 
bv iKritsuk fe Normanl (2002). We use the relation to infer 
from the observed kinematics that either the density ratio 
in the high redshift sources is different from their low red- 
shift counterparts or the temperature in the environment is 
lower. The temperature distribution function is dominated 
by a power law which we ascribe mainly to mixing. We find 
a strong peak in the distribution function at 14,000 K. This 
coincides with a strong rise in the cooling function, and cor- 
responds to an equilibrium between shock heating and adi- 
abatic compression on the one side and radiative cooling on 
the other side. We suggest that such gas is responsible for the 
optical emission in high redshifted radio galaxies, whenever 
shock ionisation dominates over photoionisation. We find a 
filling factor of 10 ~ 5 to 10 -3 for the gas in this peak. The 
velocity width for this gas increases with mass load and is 
about 100 km/s for the simulations we consider to be reliable 
in this respect. By comparison of simulations with different 
loads of cold mass, we establish that the cooling time for the 
hotter phases is considerably reduced by the presence of the 
cold gas. We find an exponential growth of the cold mass 
with time due to cooling of warmer gas. The growth rate 
is sufficient to explain the origin of the optical gas associ- 
ated with the cocoons of radio galaxies, as being build up by 
this turbulence enhanced cooling process. Heat conduction 
is not able to evaporate the cool clouds we inject. However, 
in cases where the suppression below the Spitzer value is less 
extreme than deduced from observations in nearby sources, 
it might contribute to the heat transfer significantly, and 
would presumably amplify the process we describe. 
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